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0^ ■ ABSTRACT 

I ■ A standard calculation of the energy density of heavy stable particles that may 

D I pair-annihilate into light particles making up thermal medium is performed to second 



order of coupling, using the technique of thermal field theory. At very low temper- 
atures a power law of temperature is derived for the energy density of the heavy 



> 

X. 

H , particle. This is in sharp contrast to the exponentially suppressed contribution es- 

timated from the ideal gas distribution function. The result supports a previous 
dynamical calculation based on the Hartree approximation, and implies that the 
relic abundance of dark matter particles is enhanced compared to that based on the 
Boltzmann equation. 



Estimate of the relic number density of dark matter particles such as LSP 
(Lightest Supersymmetric Particle) is conventionally made using a thermally av- 
eraged Boltzmann equation [|I|]. When the freeze-out temperature Tf is relatively 
high, for instance Tf > M/5 with M the mass of the annihilating heavy particle, 
this procedure is justified, but only after a careful examination (as we shall do in 
this paper). On the other hand, if the temperature is low, e.g, for Tf < M/30 (a 
case typical for the cold dark matter), the use of the Boltzmann equation is dubi- 
ous, and a more general quantum kinetic equation becomes necessary |0. A basic 
reason for this is that a finite time behavior of quantum system should be analyzed 
beyond the Boltzmann equation which uses the on-shell S-matrix element, hence is 
not fundamental at the full quantum mechanical level. As is well known, quantum 
mechanics at finite time cannot be described in terms of the S-matrix alone. 

Some model calculation using a kinetic equation based on the Hartree approx- 
imation was performed in 0, assuming slow variation of the occupation number. 
In cosmology the temperature variation is given by the adiabatic law, T/T = — 
the Hubble rate, and it was found that the heavy particle number density follows 
the equilibrium value until a freeze-out temperature below which the annihilation 
is frozen. This picture of the sudden freeze-out [|] is valid both for the Boltzmann 
equation and for our new kinetic equation. Thus, estimate of the freeze-out tem- 
perature using the equilibrium abundance is crucial even when one includes off-shell 
effects. It has been shown , ||^ that for the equilibrium occupation number higher 
order term in coupling dominates at low temperatures over that given by the ideal 
gas distribution function {l/{e^^'^ — 1) for bosons). A temperature power depen- 
dence was thus derived for the equilibrium heavy particle abundance. A larger relic 
abundance than previously thought of emerged, hence a more restrictive region for 
the model parameter space is anticipated. 

In the present work we employ a more familiar technique of the thermal field 
theory (imaginary-time formalism). Although there is no real time in the thermal 
field theory that governs the out-of equilibrium dynamics, the thermal average one 
computes here deals with manifestly off-shell quantities. Hence the present method 
departs from the S-matrix approach of the thermally averaged Boltzmann equation. 
(In fairness we should point out that our kinetic approach justifies the thermally 
averaged Boltzmann equation at high temperatures, but not at low temperatures.) 
We shall derive for the observable energy density a basically similar, but numerically 



different equilibrium result, from the previous result 0; the temperature power 
is different. Since definition of the occupation number is somewhat ambiguous in 
interacting field theory, a direct computation of the energy density is desirable, which 
we do in the present work. In a companion paper ^ we derive a new form of the 
kinetic equation, the one much simpler in form than that of ref. . The idea there 
is the Hartree approximation using the influence functional method, and we shall 
clarify how the detailed balance equation there agrees with the result given in the 
present work. 

For simplicity, we take throughout this paper a relativistic boson model, assuming 
the annihilation interaction 

\^'x\ (1) 

with if being the heavy particle field and x being the light field taken as making up 
a part of the thermal bath. The themal average of any operator O is done by using 
the Gibbs weight e~^^ with (3 = 1/T the inverse temperature, 

tr Oe"^^ 
(O) = -r— - (T = contribution) . (2) 

The total hamiltonian H contains both contribution from the thermal bath and 
interaction. Both for consistent renormalization and for thermalization of light x 
particles, we need self-interaction of the form, 

^^\ ^X'. (3) 

We numerically assume \X^\ to be much less than A^, but A^ of arbitrary order 
provided |A| ^ |A^| < 1. This hierarchy of coupling constants is assumed keeping 
in mind that the annihilation interaction is weak, while lighter particles can be kept 
in thermal equilibrium due to a stronger self- interaction given by A^. The mass 
hierarchy m^ <^ m^ also helps in favor of this picture. 

We compute the ip energy density p^ including renormalization counter terms; 

Pip = {Tip + (counter terms)) , (4) 

n, = -\^' + \{V^f + \M'^\ (5) 

where (p is the Euclidean time derivative. Note that by our definition of the thermal 
trace the zero-temperature value is automatically subtracted in this formula. 



To leading order 0[A°] the ip energy density given by this definition contains the 
'free field" result; 



^'^ J (27r)3 ^^/7T^/T _ 1 • ^^^ 

The quantity p^p contains product of field operators at the same spacetime point 
and require renormalization at the composite operator level, separately from the 
Lagrangian counter terms. Complication due to the operator mixing also occurs. 
For instance, the operator mixing of composite operators occurs among nine second 
rank tensor operators of mass dimension up to 4, 

Qt^v ( V5^ X^ ^\^ , ^\ X^ {9^^? , {dxf ) , <9^^<9^^ , d^xduX ■ (7) 

The renormalized Hamiltonian 7i^ is then of the form, 

^^ = ^?" + E 4'^ 02 , + E z\'^ 04 , (8) 

i i 

to 0[A^]. There are two operators 02,j of mass dimension 2, while there are seven 
operators O^^i of mass dimension 4. There is no Z^^' factor of order A due to the 
assumed interaction hierarchy. 

Renormalization constants Zi are derived by assigning a canonical value to the 
zero-temperature Fourier transformed correlation function at one external momen- 
tum point. (We use the on-shell point for renormalization.) For example, for d^ipd^^p 
the correlators with field operators up to mass dimension 4 

{o\d^ipdM^) i'i^i)i'i^2)\o) , (9) 

(Old^^dM^) ^{X1MX2MX3MX4)\0) , (10) 

with Ip either ip or x, should be taken into account. (Higher dimensional correlators 
are all finite and irrelevant to renormalization.) The free field limit values are used 
to define Z^. 

To order 0[A^] there are three types of diagrams (two topologically distinct) 
contributing to the energy density {H^). See Fig.l and Fig. 2. The temperature 
power dependence arises from the amputated sunset diagram (Fig.l), which gives in 
the configuration space 

A' / d^x dS A^{x)A^{-y)A^{y - x)A^{x - y)A^{y - x) , (11) 



prior to renormalization. The propagator in thermal medium A^^^{y) is periodic in 
the Euchdean time yo with a period (3 = 1/T, hence giving the range of integration 
~ /5 as exphcitly indicated. The Fourier transformed propagator has the well 
known form, A{uJn,p) ~ l/(— cj^ + p^ + M^) with discrete Un = 2TTin/j3 {n = 
, ±1 , ±2 ,■■■). The other contributions from Fig. 2 are exponentially suppressed by 
^-M/T ^j^]-^ ]\/[ ^]-^g heavy </? particle mass, the factor familiar in the conventional 
approach. 

As usual, one rewrites eq.(^) using the Fourier transform. The resulting discrete 
energy sum over Un can be converted to a contour integral of this variable z = Un, 
using the function l/{e^^ — 1). After some algebraic manipulation, one finds that to 

pf ~ -X^ fdk dk' dp dp' {27cf 6{p + p + k + k') 2ujp 

' UfA^ + /fc)(i + fk') - (1 + /p)(i + U')fkfk 
fpfp'ii + fk)fk' - (1 + /p)(i + U')fk{i + fl) 



+ 2 



+ 



+ 2 



+ 2 



{ujp + UJpi -UJk + UJk'Y 

fpfp'fkfk' - (1 + /p)(i + /pO(i + /fc)(i + fL) 

{iUp + up' +ujk + ^k'y 

/,(1 + /pO(1 + fk)il + f'k) - (1 + fp)fp'fkfk' 



{iOp — Up' — UJk — UJk'Y 

fpii + fp')fkii + fk') - (1 + fp)fp'ii + fk)fL 



(12) 



{Ujp ~ UJpt +UJk -UJk'Y 

We dropped minor Boltzmann suppressed terms to obtain this result. A shorthand 
notation for the phase space integral dk = d'^k/(2TT)^2uJk was used here, and fp^p/ are 
the occupation number for the heavy ip particle, while fk,k' are that for the light x 
particle; 

fp = ^,/^T^/tTI ' ^' ^ e^/^ - 1 ■ ^^^^ 

A similar form to eq.(|l2D was derived for the proper self-energy in ref.p]. For sim- 
plicity we asumme that the x mass m^ <C T, and indeed take m^ = here. 

Terms containing fp or fpi in eq.(|l^) are Boltzmann suppressed by e~^^^'^. Drop- 



ping all these Boltzmann suppressed terms, one obtains after removal of the infinity 

j^Jdk ik' U f. {^l + ^l) l^,_ is ^\/l - — . (14) 

This gives to leading order of T /M 

Terms of higher temperature-power are subleading, and neglected here. 

There is a simple reason how the temperature dependence oc T^ arises. Since 
Ti^ has both dimension 2 and 4 operators, terms of order T^ and T^ are divergent 
and they are cancelled by counter terms Zi Oi. The reason why one does not have 
divergent 0[A^T^/M^] terms is that Zi is already of order A^ and one only needs {Oi) 
to 0[A°] which is either Boltzmann-suppressed or has no T^M'"^ term. The remaining 
finite term is then of order T^ unless some special cancellation mechanism works. In 
our companion paper [Q we give a separate computation of the equilibrium energy 
density from the kinetic approach, which exactly agrees with the present result. 

We separately computed the interaction energy density given by Fig. 3 to get 

A TT^ T^ 

Pint = (^ ^\^ + (counter terms)) -^^^ A^ — , (16) 

dropping 0[A] Boltzmann suppressed terms. Thus, the interaction makes a minor 
contribution in the pair-annihilation model. 

How about terms of order AA" ? All these turn out to give Boltzmann suppressed 

A 

contributions to p^p and pint- This makes our result insensitive to the self- interaction 
among light x particles. 

A possibility that the interaction makes at low temperatures a comparable con- 
tribution to the heavy particle energy was noted by Singh and Srednicki f^ who 
explicitly calculated the interaction energy in the simple solvable model of 0. How- 
ever their suggestion that quantum kinetic approach should be discarded and we 
should go back to the Boltzmann equation is not valid for a number of reasons. First 
of all, there is no solid justification for the Boltzmann equation at low temperatures. 
Next, the amount of interaction energy is model dependent; indeed our annihila- 
tion model gives a negligible contribution of interaction. We are neither confident 
of their claim that the solvable gaussian model gives a comparable contribution of 
interaction, becasuse renormalization is not taken into account in their computation. 



Moreover, as we shall show below, one can make in the annihilation model a clear 
distinction between various forms of energy by taking into account the cosmological 
expansion. This argument makes clear which part should be regarded as the dark 
matter energy density. 

We would however like to point out some peculiarity; cancellation of order A^ T^ 
terms for the total energy, ptot = P(/3 + Px + Ant- Namely, a term in p^ of order A^ T^ 
exactly cancells the same order term in p^p to give a A^ T^, and no A^ T^, term to ptot- 
This however does not mean that both, p^ and ptot, are numerically of order A^T®. 
There are more important, larger terms of order, for instance, AA" for these. With 
a larger A^ coupling, one has a consistent picture for thermalization of x particles. 

Although it is technically difficult to compute the momentum distribution of relic 
particles, it is relatively easy to compute the pressure of heavy particles in thermal 
equilibrium. It turns out that the pressure p is one fifth of the energy, p = | p. 
It is neither of completely non-relativistic (in which case p = 0) nor of relativistic 
{p = \p) form. Thus, right after the freeze-out the equation of state implies that the 
energy density follows p oc a~ ^^'^ . But the well known redshift effect makes the high 
momentum component energetically subdominant after rapid cosmologial expansion. 
It is thus reasonable to suppose that at later epochs the heavy particle is essentially 
non-relativistic, behaving like p<^ oc a^^. On the other hand, the interaction energy 
density ~ A(/9^x^ decreases much more rapidly like a~^ in thermal equilibrium. This 
decrease is faster than the high momentum part of p^p. We thus find that the correct 
dark matter density decreasing with the volume factor should be identifed as p^. 

Since our formula for the energy density is valid only at low temperatures, T <^ 
M, meaning that the indivisual particle energy E ~ M, one can also deduce for the 
heavy particle number density at T ^ M 

Hip ^ Hip Hp . \2 ■'- t^^\ 

There is a critical temperature Tcr below which the temperature power term 
dominates over the usual Boltzmann term. This may easily be estimated by equating 
the two formulas; 

(^pe--^/^- = cA^g. (18) 

Numerically, the value of T^r/M ranges from 1/28 to 1/33 for A = 0.1 — 0.01. Thus, 
in this A range T^ ~ M/30, very crudely. A useful empirical formula of the critical 



temperature in the coupling range ofA = 0.1 — 10 ^is 

M 

— = 23 -2.3 In A. (19) 

-^ cr 

Is our result reliable at the zero temperature limit, T ^ ? We argue that this 
is not so from the following reason. Our method of using the equilibrium value for 
the freeze-out abundance is based on available sufficient time for relaxation towards 
equilibrium. In cosmology expansion makes this time limited to the Hubble time. 
The physical time scale towards equilibrium is the inverse of the pair-creation rate 

for XX -^ V^V^, 

{av)n,^0[lQ-']\\^fT. (20) 

This should be compared to the Hubble rate, 0[1] N^''^T'^ /m^i where A^ is the rel- 
ativistic degrees of freedom contributing to the energy density. The condition for 
relaxation in cosmology is then 

^ » O[10-l ^( ^ )V^ (21) 

For a very small A this condition may violate T < M, in which case the use of 
the ordinary Boltzmann equation is justified. We thus need a relatively large cou- 
pling and/or a relatively small mass M for our new result to be dominant over the 
Boltzmann suppressed number density. In our published calculation of |0] there is a 
technical mistake in computation of the if number density so that the correct con- 
dition for the off-shell effect becomes more stringent. The numerical result there 
should thus be corrected, and the parameter region for the new effect is much more 
reduced, as seen more fully in our subsequent analysis. 

We now estimate the freeze-out temperature Tj. The simplest way is to equate 
the Hubble rate given as a function of the temperature T to the annihilation rate 
{(Tv) n^p, where a is the annihilation cross section of order A^/M^. Using the number 
density, eq.(0) and equating to the Hubble rate 

H = 1MxVn , (22) 

THpi 



we find for Tj < T^ 



T,^0.3GeV^i^^r". (23) 



This holds only when the temperature power term dominates over the exponential 
^-M/T ^Q^Yo. for n^p. The relic mass density is then given by 

Demanding that this is smaller than the present critical density gives a constraint 
on the model parameter, the coupling A and the heavy mass M. The y? mass density 
relative to the critical density is at present 

^-10-^^(^^^r. (25) 

Pc A 100 Gel/ 

Numerical estimate of the freeze-out temperature including the Boltzmann sup- 
pressed region is given in Fig. 4. The relativistic degree of freedom A^ is taken 10.75 
throughout this paper. At a given A this freeze-out temperature substantially differs 
from the naive estimate obtained by using the Boltzmann factor, when the heavy 
particle mass is small. 

A more precise, yet approximate estimate is possible by using time evolution 
equation in the expanding universe. From the equation for the number density 

dn,, 



dt 



+ 3Hn^ = ~{av){n'-ntJ, (26) 



one has for the yield Y = n^/T^ 

dY 

-^ = d{av)m,,{Y'-Y^l), (27) 

\2 

d=lMVN, (av) ^ , (28) 

' ^ ^ 327rM2' ^ ' 

where the temperature-time relation t oc T^^ was used. The equilibrium value is 
approximately a sum of two terms valid at high (but T/M <^ 1/y/cX ^ 3 x 10^/A) 
and low temperatures, 

1 f d^p 1 .9.T,. 

In Fig. 5 we plot an example of numerical computation of the time evolution equation 



(|27|) . Among two examples of coupling A = 0.05 , 0.3 both with M = 100 Gel^, the 
smaller coupling case gives indistiguishable new effect and the ordinary Boltzmann 
approach is numerically correct. On the other hand, the larger coupling case gives a 
substantially different relic abundance from the conventional result. 



Time evolution obtained from eq. (^) supports the picture of sudden freeze-out 
as in the Lee- Weinberg analysis |jl[|. The major difference here is the new equilibrium 
abundance ngq. 

The final relic abundance of dark matter particles is shown as a contour plot in 
the parameter (M , A) plane. In Fig. 6 we show the present (p mass density relative 
to the critical density, for computation both with and without our new effect. Our 
new effect tends to show up for a larger coupling and a smaller mass. 

Since our new effect gives an additional positive contribution to an energy inte- 
gral, the relic density is always enhanced from the conventional result without our 
effect. Thus the allowed parameter region in the model parameter space gets always 
smaller by our result. The real question is then how much of the previously known 
region is excluded by this new effect. How our result of the relic density affects the 
presumably best motivated case of the SUSY dark matter ramains to be studied nu- 
merically. Especially in the Z and the Higgs resonance region the effective coupling 
is large [0], and there the off-shell effect may be very large. This should be checked 
in more realistic calculation beyond our boson model. 
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Figure caption 

Fig.l 

Amputated sunset diagram for the composite operator. The sohd hue and the 
broken hne are the heavy ip particle and the hght x particle propagators in thermal 
equilibrium. The crossed circle represents the insertion of the composite operator. 

Fig.2 

The rest of diagrams contributing to the composite operator as in Fig.l of 0[A^]. 

Fig.3 

Diagrams contributing to the interaction energy. 

Fig.4 

The inverse freeze-out temperature Tj^ vs the mass M of the heavy particle, 
computed for two cases of coupling, A = 0.1,0.05. Dotted curves are result using 
the Boltzmann factor. 

Fig.5 

Time evolution of the relative energy density of heavy particles, M n^/n^, com- 
puted for two cases of couphng, A = 0.3 , 0.05. For the larger couphng (A = 0.3) our 
new result and old result based on the thermally averaged Boltzmann equation do 
substantially differ, while the smaller coupling case gives indistinguishable result. 

Fig.6 

Contour plot of the present mass density of relic particles. The lines shown corre- 
spond to 1 , 10^^ , 10^^ X the critical mass density p^ (taken here 1.05x10"^ K^ GeV cm' 
with ho = 0.7), along with results obtained using the Boltzmann suppression factor. 
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